/*
 * Author: Syed Mehroz Alam
 * Email: smehrozalam@yahoo.com
 * URL: Programming Home "http://www.geocities.com/smehrozalam/" 
 * Date: 6/20/2004
 * Time: 4:46 PM
 *
 */

using System;

namespace HMMatrix
{

    public class Matrix
    {

        /*************************************************************************
         Inversion of a matrix given by its LU decomposition.

         Input parameters:
             A       -   LU decomposition of the matrix (output of LUDecomposition subroutine).
             Pivots  -   table of permutations which were made during the LU decomposition
                         (the output of LUDecomposition subroutine).
             N       -   size of matrix A.

         Output parameters:
             A       -   inverse of matrix A.
                         Array whose indexes range within [1..N, 1..N].

         Result:
             True, if the matrix is not singular.
             False, if the matrix is singular.

           -- LAPACK routine (version 3.0) --
              Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
              Courant Institute, Argonne National Lab, and Rice University
              February 29, 1992
         *************************************************************************/
        public static bool inverselu(ref double[,] a,
            ref int[] pivots,
            int n)
        {
            bool result = new bool();
            double[] work = new double[0];
            int i = 0;
//            int iws = 0;
            int j = 0;
//            int jb = 0;
//            int jj = 0;
            int jp = 0;
            int jp1 = 0;
            double v = 0;
            int i_ = 0;

            result = true;

            //
            // Quick return if possible
            //
            if (n == 0)
            {
                return result;
            }
            work = new double[n + 1];

            //
            // Form inv(U)
            //
            if (!trinverse.invtriangular(ref a, n, true, false))
            {
                result = false;
                return result;
            }

            //
            // Solve the equation inv(A)*L = inv(U) for inv(A).
            //
            for (j = n; j >= 1; j--)
            {

                //
                // Copy current column of L to WORK and replace with zeros.
                //
                for (i = j + 1; i <= n; i++)
                {
                    work[i] = a[i, j];
                    a[i, j] = 0;
                }

                //
                // Compute current column of inv(A).
                //
                if (j < n)
                {
                    jp1 = j + 1;
                    for (i = 1; i <= n; i++)
                    {
                        v = 0.0;
                        for (i_ = jp1; i_ <= n; i_++)
                        {
                            v += a[i, i_] * work[i_];
                        }
                        a[i, j] = a[i, j] - v;
                    }
                }
            }

            //
            // Apply column interchanges.
            //
            for (j = n - 1; j >= 1; j--)
            {
                jp = pivots[j];
                if (jp != j)
                {
                    for (i_ = 1; i_ <= n; i_++)
                    {
                        work[i_] = a[i_, j];
                    }
                    for (i_ = 1; i_ <= n; i_++)
                    {
                        a[i_, j] = a[i_, jp];
                    }
                    for (i_ = 1; i_ <= n; i_++)
                    {
                        a[i_, jp] = work[i_];
                    }
                }
            }
            return result;
        }


        /*************************************************************************
        Inversion of a general matrix.

        Input parameters:
            A   -   matrix. Array whose indexes range within [1..N, 1..N].
            N   -   size of matrix A.

        Output parameters:
            A   -   inverse of matrix A.
                    Array whose indexes range within [1..N, 1..N].

        Result:
            True, if the matrix is not singular.
            False, if the matrix is singular.

          -- ALGLIB --
             Copyright 2005 by Bochkanov Sergey
        *************************************************************************/
        public static bool inverse(ref double[,] a,
            int n)
        {
            bool result = new bool();
            int[] pivots = new int[0];

            lu.ludecomposition(ref a, n, n, ref pivots);
            result = inverselu(ref a, ref pivots, n);
            return result;
        }










        /// <summary>
        /// The function returns the transpose of the current matrix
        /// </summary>
        public static double[,] Transpose(double[,] a)
        {
            int rows = a.GetLength(0);
            int cols = a.GetLength(1);
            double[,] TransposeMatrix = new double[cols, rows];

            for (int i = 0; i < rows; i++)
                for (int j = 0; j < cols; j++)
                    TransposeMatrix[j, i] = a[i, j];
            return TransposeMatrix;
        }


        public static double[,] Multiply(double[,] a, double[,] b)
        {
            int a_rows = a.GetLength(0);
            int a_cols = a.GetLength(1);
            int b_rows = b.GetLength(0);
            int b_cols = b.GetLength(1);

            if (b_rows != a_cols)
                throw new Exception("Matrix: Operation not possible");

            double[,] result = new double[a_rows, b_cols];
            for (int i = 0; i < a_rows; i++)
                for (int j = 0; j < b_cols; j++)
                    for (int k = 0; k < a_cols; k++)
                        result[i, j] += a[i, k] * b[k, j];

            return result;
        }

    }//end class Matrix

}	//end namespace HMMatrix
